Escape transition of a polymer chain from a nanotube: how to avoid spurious results 
by use of the force-biased pruned-enriched Rosenbluth algorithm 



Hsiao-Ping Hsu and Kurt Binder 

Institut fur Physik, Johannes Gutenberg- Universitat Mainz 
D-55099 Mainz, Staudinger Weg 7, Germany 



Leonid I. Klushin 

American University of Beirut, Department of Physics, Beirut, Lebanon 

Alexander M. Skvortsov 
Chemical-Pharmaceutical Academy, Prof. Popova 14, 197022 St. Petersburg, Russia. 

(Dated: August 11, 2008) 

A polymer chain containing N monomers confined in a finite cylindrical tube of diameter D grafted 
at a distance L from the open end of the tube may undergo a rather abrupt transition, where part of 
the chain escapes from the tube to form a "crown-like" coil outside of the tube. When this problem 
is studied by Monte Carlo simulation of self-avoiding walks on the simple cubic lattice applying 
a cylindrical confinement and using the standard pruned-enriched Rosenbluth method (PERM), 
one obtains spurious results, however: with increasing chain length the transition gets weaker and 
weaker, due to insufficient sampling of the "escaped" states, as a detailed analysis shows. In order 
to solve this problem, a new variant of a biased sequential sampling algorithm with re-sampling is 
proposed, force-biased PERM: the difficulty of sampling both phases in the region of the first order 
transition with the correct weights is treated by applying a force at the free end pulling it out of the 
tube. Different strengths of this force need to be used and reweighting techniques are applied. Using 
rather long chains (up to TV = 18000) and wide tubes (up to D — 29 lattice spacings), the free energy 
of the chain, its end-to-end distance, the number of "imprisoned" monomers can be estimated, as 
well as the order parameter and its distribution. It is suggested that this new algorithm should be 
useful for other problems involving state changes of polymers, where the different states belong to 
rather disjunct "valleys" in the phase space of the system. 



I. INTRODUCTION 



A polymer chain with one end grafted to a non- 
adsorbing impenetrable surface exists in a mushroom 
conformation. If it is progressively squeezed by a flat pis- 
ton of finite radius L from above the conformation gradu- 
ally changes into a relatively thick pancake. However, be- 
yond a certain critical compression, the chain configura- 
tion changes abruptly. One part of the chain forms a stem 
stretching from the grafting point to the piston edge, 
while the rest of the segments forms a coiled crown out- 
side the piston, thus escaping from the region underneath 
the piston. This phenomenon was named the escape tran- 
sition and has attracted great interest [j], 0, HI, U HI • An 
abrupt change from one state to another implies a first 
order transition. Phase transitions at the level of an in- 
dividual macromolcculc have generic features that are 
common to all conventional phase transitions in fluids 
or magnetics. On the other hand, they may be quite 
unusual, and the conceptual framework that would en- 
compass all the specifics of this class of phase transitions 
is still being elaborated. A single macromolecule always 
consists of a finite number of monomers N: computer 
modeling rarely deals with N larger than 10 5 [6| so that 
finite-size effects in the single-molecule phase transitions 
are the rule rather than exception. The concept of a 
phase transition requires that the thermodynamic limit 
is explored. For a single macromolecule this means the 



limit of Na — > 00, while in the specific case of the escape 
transition it was shown that the appropriate limits are 
Na — > 00 and L — > 00 but Na/L = const. 

The classical examples of phase transitions in a single 
macromolecule are the coil-globule and adsorption transi- 
tions as well as the coil-stretch transition in a longitudinal 
flow that have been studied for many years theoretically, 
experimentally and by computer simulations. A dramatic 
progress in experimental methods allowing detection and 
manipulation of individual macromolecules happened in 
the last decade and includes the AFM manipulations, op- 
tical tweezers, high- resolution fluorescent probes [3, [a Si • 
As for the escape transition in its standard setup, ensur- 
ing that the piston is flat and parallel to the grafting sur- 
face currently turned out to be a difficult technical prob- 
lem. The progress in these methods is, however, so rapid 
and impressive that experimental studies of the escape 
transition must soon become within reach. Meanwhile 
the theoretical analysis suggests that the escape transi- 
tion is extremely unusual in many aspects as compared 
to conventional phase transitions in fluids or magnetics. 
First, this is a first-order transition with no phase coexis- 
tence, and therefore, no phase boundaries, no nucleation 
effects, etc. On the other hand, metastable states are 
well defined up to spinodal lines; they have a clear phys- 
ical meaning and can be easily visualized. The transi- 
tion is purely entropy-driven, the energy playing no role 
at all. Of course, this is a consequence of considering 
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FIG. 1: Schematic drawings of a flexible polymer chain of length N grafted to the bottom of the left side of which is open at 
the right end, of a tube of length L and diameter D at the transition point. As the chain is fully confined in the tube (in an 
imprisoned state), it forms a sequence of rib = ND -1 ^ blobs in a cigar-like shape, here v is the 3D Flory exponent. As one 
part of the chain escapes from the tube (in an escaped state) , it forms a flower-like configuration which consists of a "stem" 
and a coiled "crown". By adjusting one of the three parameters D, L and N, an escape transition occurs. At the transition 
point D = D* , L — L* , or N — N* , the conformation of a polymer chain is either in an imprisoned state or in an escaped 
state. 



only the excluded volume and hard-wall potentials, and 
the similar entropy-driven transitions were found in the 
hard-sphere and hard-rod fluids [lfj. Further, the escape 
transition demonstrates negative compressibility and a 
van-der-Waals-type loop in a fully equilibrium isotherm: 
this follows from exact analytical theory and was verified 
by equilibrium Monte-Carlo simulations. Finally, it was 
shown that the equation of state (compression force ver- 
sus the piston separation) corresponding to the escape 
transition setting is not the same in two conjugate en- 
sembles (constant force and constant separation) giving 
a unique example of the non-equivalence of statistical en- 
sembles that persists in the thermodynamic limit. Natu- 
rally the literature on the subject is quite extensive 

Yet one more intriguing feature of the escape transi- 
tion is that it persists in quasi-one-dimensional settings 
such as a 2d chain confined in a strip or a 3d chain con- 
fined in a tube. The existence of such a transition seems 
to formally contradict a well-known statement forbidding 
phase transitions in Id systems [l2[ . It is also counter- 
intuitive and goe s ag ainst a simple blob picture of a chain 
in a tube [11,13) Elllj[i3> which suggests that squeez- 
ing the tube would result in progressively pushing out 
the chain segments very much alike the toothpaste. Pre- 
dictions of the naive blob picture were discussed in the 
context of a 2d chain confined in an open strip of finite 
length 18]. It was shown that they contradict both the 
MC simulations with PERM and a more sophisticated 
analytical theory. There are two reasons that prompt 
a revisiting of the escape transition in one-dimensional 
geometry. First, the strip confinement of a 2d chain is 
rather exotic and difficult to realize experimentally. On 
the other hand, confining a real chain in a nanotube with 



the diameter much smaller than the chain gyration radius 
in solution is experimentally feasible. Well-calibrated 
nanochannels are being produced by lithographical meth- 
ods with the width in the range between 30 and 400 
nm (l9j . Partially escaped configurations appear also 
in situations where a long polymer chain is translocat- 
ing through a pore in a relatively thick membrane [20l |. 
From the point of view of computational physics, in gen- 
eral, first-order transitions and inhomogeneous polymer 
conformations are particularly very difficult to cope with; 
we have encountered these difficulties in the study of a 2d 
escape transition [l8j], which motivated the development 
of a new variant of the algorithm PERM, force-biased 
PERM. This new algorithm is specifically suited for ex- 
ploring well-separated regions in the phase space. 

The aim of our paper is twofold: (1) we demonstrate 
unambiguously the existence of a first-order escape tran- 
sition in the setting of a real 3d polymer confined in a 
tube; (2) we present the force-biased PERM, and show 
its advantages in comparison with PERM, or more pre- 
cisely, PERM with k-step Markovian anticipation. 

The paper is organized in a non-traditional way. In 
Sec. II, we present three different scenarios of the escape 
transition and the simulational model. In Sec. Ill and 
Sec. IV, we first present the results of the Monte Carlo 
(MC) simulations with PERM and force-biased PERM. 
These include a collection of equilibrium characteristics 
and lateral size distribution for a chain in an infinite tube 
(without the possibility of escape) in comparison with the 
scaling predictions. Then we discuss the characteristics 
of the escape transition and its smoothing due to the 
finite-size effects. In Sec. V. and VI, we explain the idea 
and the implementation of force-biased PERM, and give 
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FIG. 2: Schematic drawing of a one-end grafted polymer chain growing as a self-avoiding walk inside a finite tube of length L 
and diameter D along the x-direction. The first monomer is located at (x — 0, y = 0, z = 0). Other monomers are allowed to sit 
on lattice sites of a simple cubic lattice, except for the lattice sites representing the cylindrical walls {0 < x < L, y 2 + z 2 = D 2 /4} 
and the bottom {a; = 0, y' 2 + z 2 < D 2 /4}. 



a clear demonstration of a failure of PERM in applica- 
tion to the escape transition. We also present a glimpse of 
the computational technique whereby the artifacts of the 
method could lead to misleading conclusions concerning 
the nature of the transition in the thermodynamic limit, 
and discuss ways to detect the possible flaws in the sim- 
ulation results. Finally, a summary is given in Sec. VII. 

II. ESCAPE TRANSITION SCENARIOS AND 
SIMULATIONAL MODEL 

The escape transition of an end-attached polymer 
chain inside a confined space may be triggered by chang- 
ing any of the following parameters: 

(a) the piston separation [H,I2> 0,0,1^] or the tube diam- 
eter [ijj which control the degree of confinement; 

(b) the distance from the grafting point inside the con- 
fined space to the edge of the opening (this would 
be the piston radius R in a traditional escape set- 
ting [H, [I, d, H [H or the tube length L) 

(c) the chain contour length, Na, where a is the seg- 
ment length, a is set to 1 throughout the paper 
hereafter. 

Three corresponding scenarios of a polymer chain es- 
caping from a tube are schematically illustrated in Fig. [I] 
although not all of them can find a reasonable exper- 
imental implementation. In a classical escape setting, 
the confinement width is easily changed by moving the 
pistonjwhile the piston radius was assumed to be a con- 
stant P, 0, IH • For a tube, changing its diameter 
D is much more difficult to realize experimentally, but 
the distance from the grafting point to the tube opening, 
L, can be changed naturally by moving a bead with a 
chain end grafted to it, inside the tube. Chopping the 
tube off in thin slices is another way of changing L albeit 
rather artificial. Finally, one can envisage a mechanism 



which leads to a gradual change in the chain contour 
length (this may involve ferments that cut a molecule). 
From the point of view of PERM, the scenario of gradual 
growing a chain is the easiest one to be implemented. We 
study single flexible polymer chains of length N with one 
end grafted to the center of the bottom of a finite tube 
of length L and diameter D under good solvent condi- 
tions. They are described by a self-avoiding random walk 
(SAW) of N steps on a simple cubic lattice as shown in 
Fig. [2] with the restriction that monomers are forbidden 
to be located on the surface of the cylinder. Taking the 
bottom center at the origin (x, y, z) — (0, 0, 0) with the 
first monomer fixed there, and the tube axis along the 
x-direction, monomers are not allowed to be located on 
the lattice sites < x < L, y 2 + z 2 = D 2 /4 and x = 0, 
y 2 + z 2 < D 2 j\. The idea of using a half-opened tube 
and a chain with one end grafted to its bottom is that 
with a chain growth algorithm the chain has to grow ei- 
ther along +x-direction or — ^-direction as N v ~ D, and 
the partition sum just differs by a constant comparing to 
the case of unfixed chain ends [2lJ. However, the scaling 
behaviour of confined polymer chains remains the same 
for this simplified model. 



III. CONFINED CHAINS IN AN INFINITE 
TUBE 

The properties of a single macromolecule confined 
in a tube have been studied extensively for decades, 
both by analytical theory and by numerical simula- 
tions for various models of flexible and semi-flexible 
chains [H, iii For a homogeneous confined 
state there are scaling predictions fl3j concerning various 
chain characteristics which were tested by MC simula- 
tions. 

The scaling laws are well known [13 . |22| for a fully 
confined polymer chain of chain length TV in a tube of 
length L and diameter D under a good solvent condition. 
One should expect a crossover between a weak confine- 
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ment region 1 <C <C D (y is the Flory exponent), 
where the chain forms a three-dimensional random coil, 
and a strong confinement region 1 <C D <C N", where 
the chain forms a quasi-one dimensional conformation (a 
long cigar-shaped object which is described as a chain of 
spherical blobs in a blob picture). 

For a grafted polymer chain confined in a very wide 
tube, D » N v , since the chain stays essentially unaf- 
fected by the existence of the cylindrical hard surface, 
i.e., there is no repulsive interaction between the chain 
and the cylindrical hard surface, the system is simplified 
to be in the situation that a polymer chain is with one 
end grafted to a repulsive wall in a good solvent. Using 
a self-avoiding walk (SAW) of N steps on a simple cubic 
lattice with a constraint that x > 0, it is well known that 
in the limit of N — > oo the end-to-end distance scales as 



R N ~ N" 



(1) 



with v = 0.58765(20) [23|], which shows the same behav- 
ior as in the bulk. 

A scaling ansatz of the end-to-end distance along the 
direction parallel to the tube axis, R\\(N, D), for the 
cross-over between these two regimes 1 << N v « D 
and 1 << D « N v is 



Ru(N,D)=R N $ R (R N /D) 



(2) 



where Rn is given by Eq. JT]) and requesting that for 
large N the parallel linear dimension is proportional to 
N yields for the function (f>^ the limiting behaviours 
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(3) 



with r] — N v /D. The same scaling ansatz is also ex- 
pected for R X (N,D), 



R ± {N,D) = R N V R {R N /D) 



(4) 



and requesting that for large N the perpendicular linear 
dimension is independent of N yields for the function ^ r 
the limiting behaviours 




N v /D 

FIG. 3: Rescaled rms end-to-end distance parallel to the tube, 
R\\(N,D)/N", plotted against n = N v /D. The data collapse 
is achieved by Eq. ((2} and corrections to scaling are not vis- 
ible in a log-log scale. The inset shows the log-log plot of 
A-Riioc = linijv^oo R\\(N,D)/N versus D. The dashed line is 
0.93D 1_1/l/ . 
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FIG. 4: Rescaled rms end-to-end distance perpendicular to 
the tube, R± (N, D) /TV" , plotted against rj. The data collapse 
is achieved by Eq. (j4]) and corrections to scaling are not visible 
in a log- log scale. The inset shows the log- log plot of R±,oo = 
limjv^oo R ± (N,D) versus D. The dashed line is 0.295D. 
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(5) 



In order to obtain the full functions $7^(77), ^air/) 
interpolating between the quoted limits, we apply MC 
methods. For the MC simulations, we use the chain 
growth algorithm PERM with fc-step Markovian antic- 
ipation and simulate single fully confined chains of chain 
lengths N up to 44000, and tube diameters D up to 97. 
According to the cross-over scaling ansatz, Eqs. {2} and 
(jl]) , we plot the rescaled average root mean square (rms) 
end-to-end distances parallel and perpendicular to the 
tube axis, i?||/A l/ and R±/N u , against N"/D in Fig. [3] 
and HI respectively. A nice data collapse is indeed seen 
in both figures. For a further check of the scaling law 
depending on D in the regime of N v » D (the strong 



confinement limit) and giving a precise estimate of the 
prefactor, we estimate the asymptotic ratios between 
R\\ and A, i.e. AR\\ tOC = ]hn N ^ 00 R\\{N,D)/N. and 
the asymptotic value of R± as N — > 00, i.e. R±.oo — 
limTv^oo R±{N, D), for various diameters D. Results of 
Ai2n i0O and Rx,oot plotted against D are shown in the 
inset of Fig. [3] and Fig. [H respectively. In the thermody- 
namic limit N — > 00, Ru and Ri scale as follows 



Ru — 



0.85(2)A iy 
Q.93(2)ND 1 - 



for D » N u » 1 
V" for N v »D»1 



and 



R± = 



0.89(2) N u for D » N v » 1 
0.295(8)L» for N v » D»l 



(6) 



(7) 
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Due to the confinement of the chain in a tube, in Eq. ([7]) 
we see that R± is only dependent on the tube diameter D 
but independent of N in the strong confinement regime. 
In the very weak or no confinement regime, the estimate 
of the ratio between the average mean square end-to-end 
distances in both directions, (R\\/R±) 2 m 0.91(8) is in 
good agreement with the result in Ref. [2J] (< z 2 N > / < 
x n + Vn >= 0.938(2) as TV — > oo), although our chain 
lengths in the region D >> N v >> 1 are much shorter. 

It is also remarkable that the crossover between both 
limits in Eqs. ©, Q is rather sharp and not spread out 
over a wide regime in the scaling variable N u / D. For R\\ 
this crossover occurs for N u / D w 1 but for R± it occurs 
near N v /D w 0.4, however. 

An additional bonus of the PERM algorithm 
as compared to standard "dynamic" MC algorithms 
as applied in [3, HH is that it yields the free en- 
ergy. The excess free energy relative to a one-end 

grafted random coil, F(N, D) - 



In 



Z(N,D)/Z$ 



where Z\ 



,-N 



- H *'iVT 1 with fi = 0.21349098(5) 
and 7W = 0.6786(12) is a constant in the regime 
D >> N u >> 1 while it is proportional to the number 
of blobs, nb — ND -1 ^, in the unit of fc^T in the regime 
N u >> D >> 1. Therefore, we can rewrite the scaling 
law as follows, 



F(N, D) _ f const n' 1 ^ for n -> 
ND^ 1 /" \ const' for n — > oo 



(8) 



As expected, results of the rescaled free energy in Fig. [5] 
also show the nice data collapse of the cross-over behavior 
predicted by Eq. ([8]). In the strong confinement regime, 
the imprisoned free energy scaled as follows, 



F imp (N,D) = 5.40(2)n 6 



(9) 



Again the crossover is very sharp, occurring near a value 
N u / D w 0.4 of the scaling variable. 

According to the scaling theory the free energy 
in the strong confinement limit is always proportional to 
the number of blobs, n^. In the special case of quasi-ld 
confinement (a tube for 3d chains or a strip for 2d chains) 
the average end-to-end distance is proportional to n^D, 
so that the ratio FD/R\\ is independent of both N and 
D. Physically, this ratio gives an estimate for the free 
energy of confinement per blob in ksT units. The ques- 
tion of whether this ratio is model-dependent or universal 
was addressed by Burkhardt and Guim [25| . Using field- 
theoretical methods and the equivalence of self-avoiding 
walks and the n- vector model of magnetism in the limit of 
n — > they have shown that this ratio (which we denote 
in the following as the Burkhardt amplitude Ab) is in- 
deed universal. Extrapolating finite-size transfer-matrix 
results for the case of repulsive walls in d = 2 they ob- 
tained a numerical estimate of Abi ~ 2.10 ± 0.01. Al- 
though MC simulation data for 2d SAWs confined in a 
strip with hard walls were obtained a few years ago, no 
direct comparison to the theoretical prediction was made 
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FIG. 5: Rescaled free energy relative to a one-end grafted 
random coil, F(N, D)/{ND~ 1/ ") in a log-log scale. The data 
collapse is achieved by Eq. Q. The inset shows the log- 
log plot of AF im p,oo = limiv-K» F{N,D)/N versus D. The 
dashed line is 5A0D' 1/v . 



at the time. The MC data of Hsu and Grassberger give 
i?H w 0.915AD" 1 / 3 [13 and F = 1.944(2)A^Z)- 4 / 3 |l| 
which results in Ab2 ~ 2.12(4) in excellent agreement 
with field-theoretical calculations. For a 3c? chain in a 
tube the Burkhardt amplitude must be also universal al- 
though no numerical estimates were ever produced. From 
our MC results, Eq. © and ©, it gives A B3 » 5.79(6). 
It would be interesting to check this value using other 
simulation models. 

It was proposed and checked numerically for short 
polymer chains in Ref. [l6| that the distribution for the 
gyration radius r g along the tube axis can be presented 
as a sum of two terms: 

\nP{r g \N,D) = N{D/ a y 1/u 'A[ u - a + Bu s ] ,(10) 

where a and 5 are linked to the space dimension d and the 
Flory exponent v by a = [yd — 1) 1 and S = (1 — v) 1 , 
andu = (r g /Na){D/a)- 1+1 / v is the segment volume con- 
centration expressed as a function of the gyration radius 
and the confinement geometry. The parameters A and 
B are model-dependent numerical coefficients of order 
unity, which do not depend on N or D. 

From our MC simulations, we found that the above for- 
mula has to be corrected by an additional r g -independent 
term, i.e. 



lnP(r g \N,D) = N(Dfa)^ 1 ^ A [u~ a + Bu s + C] . 

(11) 

Results of the rescaled distribution 
ln(P(r g \N,D)/P{r m \N,D))/ND- 1 ^ plotted against 
r g /(ND 1 ~ l / v ) for N = 5000, and 10000, and for 
D = 25, 33, and 49 are shown in Fig [6] Here the 
radius of gyration r g ^ m corresponds to the maximum of 
the distribution P(r g \N,D). Near r g ^ m , we see a very 
nice data collapse, and the rescaled distribution can be 
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FIG. 6: Rescaled distribution 

ln(P(r g \N,D)/P(r g , m \N,D))/(ND- 1/u ) plotted against 
rg/ND 1 - 1 '" fo r various values of N and D. The radius of 
gyration r g , m corresponds to the maximum of the distribu- 
tion P(r g \N, D). Near r g , m , f(x) = A(x~ a + Bx s + C) with 
A = -0.185, B = 76.757, and C = -8.784 gives the best fit 
of the data, here x = r g /{ND 1 - 1 / v ). 




r,, / (ND 1 " 1/V ) 



FIG. 7: Rescaled distribution 

\n{P(r ll \N,D)/P{r llrn \N,D))/{ND~ 1/,/ ) plotted against 
ru/ND 1-1 '" for various values of N and D. The end-to-end 
distance r« m corresponds to the maximum of the distribution 
P(r\\\N,D). Near r g , m , f(x) = A(x- a + Bx s + C) with 
A = -0.953, B = 0.801, and C = -1.771 gives the best fit of 
the data, here x = r y / '(N 'D l ~ 1/v ). 



described by Eq. dTTJ with A = -0.185, B = 76.757, and 
C = -8.784 very well. Thus we conclude that Eq. (fTTTl . 
even if it is not exact, at least is a very good numerical 
approximation to the actual distribution function. In 
Fig. we check that the results for the end-to-end 
distance distribution P(r\\\N, D) can be also described 
by the same function {Eq. 1111} with A = —0.953, 
B = 0.801, and C = -1.771. 




ND" 1/V 



FIG. 8: Free energy relative to a one-end grafted chain, 
F(N, L, D) = -]n[Z(N,L,D)/Z$ ) (D)], plotted against 
ND" 1/v . 



IV. SINGLE POLYMER CHAINS ESCAPE 
FROM A NANOTUBE 

We simulate 3D SAW and BSAWs starting at the graft- 
ing point of the tube with length L = 200, 400, 800, and 
1600, and diameter D = 17, 21, 25 and 29. Depending 
on the chosen size of L and D, the total chain length is 
varied from 1400 to 18000 in order to cover the transition 
regime. Results of the free energy relative to a one-end 
grafted random coil for L — 800 and L — 1600 are shown 
in Fig. [51 It is clear that there are two branches of the 
free energy. For a fixed diameter D we see that initially 
the free energy increases linearly with the chain length 
N, but as N exceeds a certain value, the chain escapes 
from the tube, a sharp crossover behavior from an im- 
prisoned state to an escaped state is indeed seen. Values 
of the excess free energy of the chain in an escaped state, 
F esc (N, L, D), are determined by the horizontal curves 
shown in Fig. [8j which are independent of N. Results of 
F esc (N, L, D) for various values of L and D are presented 
in Fig. [SI where we obtain 

F esc (N,L,D)=i.23(7)L/D (12) 

According to the requirement that the free energies 
of polymer chains in both imprisoned and escaped states 
should be the same at the transition point, the transition 
point is therefore determined by equating Eq. ([9]) and 
Eq. (fT2"]) . We obtain the relation between L, N, and D 
at the transition point, i.e., 

(N/L) tI = 0.78(2)D 1 / y - 1 . (13) 

One can also estimate ND -1 ^ at the transition point 
directly from Fig. [SJ for fixed tube length L and tube 
diameter D. In Fig. [10] we plot ND~ 1/u against L/D. 
The best fit for the data given by the straight line is 
ND~ X I V = 0.785(10)L/-D which is in perfect agreement 
with Eq. ([13}. 
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FIG. 9: The log-log plot of the excess free energy of the es- 
caped chain, F eac (N, L, D), plotted against L/D for various 
values of L and D. The dashed line is 4.23L/D and gives the 
best fit of the data. 



FIG. 11: Average end-to-end distance divided by L, R\\/L, 
plotted against N/L for various values of L and D. The 
rounding of the transition is due to the finiteness of the chain 
length N. 
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FIG. 10: Results of the transition points obtained by esti- 
mating the intersection points of the two branches shown in 
Fig. IS] The dashed line is 0.785L/D and gives the best fit of 
the data. 




FIG. 12: Average fraction of imprisoned number, Ni mp /N, 
plotted against N/L for various values of L and D. The 
rounding of the transition is due to the finiteness of the chain 
length N. 



In order to understand the conformational change of 
the polymer chains from an imprisoned state to an es- 
caped state, in Fig. [11] we show the results of the end-to- 
end distance parallel to the tube axis normalized by the 
tube length L, Ru/L, versus N/L. As the chain is still in 
the strong confinement regime, Ru increases linearly with 
N, showing that the chain is stretched. Beyond a certain 
value of N/L, there is a jump in each curve, showing the 
chain undergoes an escape transition, i.e., one part of 
the chain escapes from the tube (see Fig. [TJ. Obviously, 
we see that the transition becomes sharper as the tube 
length L increases or the tube diameter D decreases. The 
pronounced jump indicates that the transition is first- 
order like. A similar jumpwisc change is also expected 
for the average number of imprisoned monomers, N mp . 
In Fig. 1121 we plot the fraction of imprisoned monomers, 
Nimp/N, versus N/L for various values of L and D. It 
is indeed seen that there exist jumps from N mp /N = 1 



while the chain is confined completely to N mp /N * 3/4 
as N/L increases, and the rounding of the transition is 
due to the finite size effect. In the problem of escape 
transition, the order parameter is defined by the stretch- 
ing degree of the confined segments of polymer chains. 
As the chain is in an imprisoned state, s = r/N where r 
is the instantaneous end-to-end distance of the confined 
chain along the tube axis, while as the chain is in an 
escaped state, s = L/n where n is the imprisoned num- 
ber of monomers (number of monomers in the stem) [26| . 
Results of the average order parameter S =< s > against 
N/L are presented in Fig.[T2J An abrupt jump from S- lmp 
to S esc is developed at the transition point in Fig. [T2] as 
the system size L increases for each data set of a given 
value of D. Since it is difficult to give a precise estimate 
of the transition point directly from those results shown 
in Fig. [TTIfTSI we use a different strategy by plotting a 
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FIG. 13: Average order parameter S plotted against N/L for 
various values of L and D. The rounding of the transition is 
due to the finiteness of the chain length N. 
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FIG. 14: Distribution functions of the order parameter 
P(N, L, D, s) vs. s near the transition point for L = 1600 
and D = 17. . The distribution is normalized by choosing 
^ s P(N, L, D, s) — 1. A constant A is added to distinguish 
between different curves. 



straight line with N/L = Q.785D 1 /"- 1 for each D. We 
see that the straight lines do go through the intersection 
point for a fixed value D (results are not shown) . 

Finally, results of the distribution function of the or- 
der parameter, P(N 7 L,D,s), are presented in Fig. [T4l 
for L — 1600 and D — 17 near the transition point. For 
our simulations, the distribution function P(N, L, D, s) 
is given by properly normalized accumulated histograms 
of the order parameter s, namely J2 S L, D, s) = 1 
(see Sec. V). We see the bimodal behavior of the distri- 
bution functions as one should expect for the first-order 
transition near the transition point. At the transition 
point (N/L) tT w 5.7, the two peaks are of equal height. 
Below the transition point chains are in favor of staying 
in an imprisoned state, while above the transition point 
chains stay in an escaped state. 



FORCE-BIASED PERM 



According to the predictions of the Landau free en- 
ergy approach in three dimensions at the escape tran- 
sition point, the magnitude of the jumps of the end-to- 
end distance, the order parameter, and the fraction of 
imprisoned number of monomers should be much more 
prominent than that for the 2d case described in our 
previous work [j^ . There, polymer chains were de- 
scribed by 2d SAW on a square lattice with the con- 
straint that chains are confined in finite strips. In order 
to suppress dense configurations and sample more rela- 
tively open chain configurations, we employed the algo- 
rithm PERM with fc-step Markovian anticipation [H, [2l[ 
to study the problem of 2d escape transition. However, 
this method is very efficient for generating homogeneous 
configurations of very long polymer chains in the impris- 
oned state not only for the 2d case but also for the 3d 
case which will be presented in the next section, while it 
fails with producing inhomogeneous flower-like configu- 
rations in the escaped state. The difficulty was already 
noticed in Ref. [18] by the disappearance of jumps in the 
thermodynamic limit and the lack of data for describing 
the Landau free energy as a function of the order param- 
eter s in an escaped state for larger systems. Clearly, it 
becomes a more serious problem when the same sampling 
method and model are applied to the study of polymer 
chains escaping from a tube. Our test run showed that for 
a fixed diameter D the transition occurs much later (if we 
choose N/L as a control parameter) and the magnitude 
of the jump decreases as the tube length L increases (see 
Fig. I15p . One could imagine that in the thermodynamic 
limit of N — > oo and L — > oo, the conformation change 
becomes continuous. However this conjecture is based on 
the artifacts due to the inefficiency of the original chosen 
model and method. 

A new strategy for generating sufficient samplings of 
the flower-like configurations in the phase space is pro- 
posed as follows: We first apply an extra constant force 
along the tube to pull the free end of a grafted chain out- 
ward to the open end of the tube as long as the chain 
is still confined in a tube, and release the chain once 
some segments of it is outside the tube. By varying the 
strength of the force, we obtain configurations in the es- 
caped state with various stretching degree of monomer 
segments (stems) which are still confined in a tube. Fi- 
nally the contributions for the escaped states are given 
by properly reweighting these configurations to the case 
without applying extra forces. Now a partially stretched 
polymer in a good solvent is described by a biased SAW 
(BSAW) with finite cylindrical geometry confinement. 
The stretching is denoted by a factor b — exp((3aF) where 
a is the lattice constant, F is the stretching force, and 
(3 = 1/ksT. For our simulations, a and (3 are rescaled 
as units of length and inverse energy, respectively. Both 
are set to 1 in the simulations. The partition sum of a 
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BSAW of N steps is, therefore, 

Z b (N,L,D) = Y, b 



Ax 



with 



> 1 , < x < L ,y 2 + z 2 < D 
1 , otherwise (SAW) 



(14) 



(15) 



here Ax = xm — xq is the displacement (in units of a lat- 
tice constant) of the end-to-end vector onto the direction 
of F (along the tube axis). The first monomer is located 
at xq = as shown in Fig. [5J Based on the algorithm 
PERM, polymer chains are built like random walks by 
adding one monomer at each step and each configura- 
tion carries its own weight which is a product of those 
weight gains at each step, i.e. W(N) = n|^ 1 w(j) with 

w (j) = b^ - *^ -1 ) and w(j = 0) = 1 in the current case 
(cf. Eq. (HU). It has the advantage that one can estimate 
the partition sum directly as given by, 



Z b (N,L,D) = l- J2w b (C b ) 



(16) 



{Cb} 



here C b denotes a configuration of a BSAW of N steps, 
confined in a finite tube of length L and diameter 
D, M b is the total trial configurations, and W b (C b ) — 
W b (N, L,D) is the total weight of obtaining the configu- 
ration C b . Thus, each configuration of a BSAW with the 
stretching factor contributes a weig ht W {k \N, L,D) 
after re-weighting to compensate for the bias introduced 
earlier as 

W^(N L D)-( W^N.^D)^-^ , x N < L 
[ ' ' U) -\W bk (N,L,D)/bl , x N >L ' 

(17) 

here index k labels runs with different values of the 
stretching factor b. 

Combining data runs with different values of b the av- 
erage value of any observable O is given by 

,^_ ^ C on flg .ec bk O(C bk )W^(N,L,D) 

E k J2 conflg .ec bk W^(N,L,D) [1 *> 

and the estimate of the partition sum 

Z(N,L,D) = ±-Y1 E W™(N,L,D) (19) 



k config.^Cb,, 



here M is the total number of trial configurations. 

The distribution of the order parameter, 
P(N,L,D,s) oc H(N,L,D,s), is obtained by accu- 
mulating the histograms H(N,L,D,s) of s, where 
H(N, L, D, s) is given by, 

H (N,L,D,s) 

k 

JjJ2 E WW(N,L,D,s')5 StS , (20) 



k config.£Cb k 



and the partition sum of polymer chains confined in a 
finite tube can be written as 

Z(N,L,D) = J2H{N,L,D,s) (21) 

s 

in accordance with Eq. (|19p . 



VI. COMPARISON OF RESULTS BY OLD AND 
NEW ALGORITHM 

Here we are going to show some of the technical de- 
tails of simulations by PERM and the resulting problem. 
We present the spurious results obtained by PERM with- 
out force biases in detail and the far-reaching conclusions 
they could lead to. The whole situation is of considerable 
methodological and pedagogical value. We also discuss 
some general ways to avoid the pitfalls and discriminate 
between authentic and spurious data. 



A. Results without force biases: new phase 
transition physics suggested? 

Let us focus on the results of the average end-to-end 
distance R\\/L. Using the algorithm PERM with fc-step 
Markovian anticipation as in our previous work (2d es- 
cape transition) [13], we present the data of R\\/L as a 
function of the reduced chain length parameter N/L for 
D = 21 in Fig.[TSJ Each curve corresponds to a fixed tube 
length L. With increasing the system size, i.e. increas- 
ing L from L = 200 to L = 1600, the systematic change 
of the curves describes the approach to the thermody- 
namic limit. The data points fall on smooth curves with- 
out much statistical scattering, and the curves themselves 
suggest a systematic trend as L increases. It is clear that 
for each L the chain size experiences a sharp increase at 
some value of N/L. The sharpness of the curves increases 
to some extent as the system size increases near the tran- 
sition point, which is not unexpected if one suspects a 
first-order transition to be involved. The transition point 
of a finite system can be estimated, and a shift of its 
position with increase in L is obvious. Thus one is im- 
mediately tempted to extrapolate the transition point to 
L — > co, which is pointed by an arrow in Fig. 1 151 Simulta- 
neously, the magnitude of the jump is also changing sys- 
tematically which calls for another extrapolation. This 
extrapolation shows that in the thermodynamic limit the 
jump vanishes (the long dashed line in Fig. 1151 ) Overall, 
this suggests a sophisticated and weak first-order transi- 
tion with finite-size effects that disappears in the ther- 
modynamic limit. This is in agreement with the intu- 
itive picture of a chain overgrowing and escaping out of 
the tube opening without any jumps and also follows 
the conventional blob picture. Further analysis would 
have to reconcile several conflicting findings. At the ex- 
trapolated transition point, the free energy still has a 
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FIG. 15: Average end-to-end distance divided by L, R^/L, 
plotted against N/L for various values of L and D — 21. 
Results are obtained by PERM without force biases. The 
transition point indicated by the arrow is the intersection be- 
tween the extrapolated curve (long dashed curve) and the 
curve R\\/L = 1. 



ases for L = 200 and 1600, for D — 21, and for the chain 
length N which is large enough to cover the transition 
regions are shown in Fig. \W\ by symbols and curves re- 
spectively for comparison. It is clear that for the smaller 
system, L = 200, both algorithms give the same results. 
In contrast to this for L = 1600 the new algorithm pro- 
duces a very different curve that does not show a shift 
of the transition point (nor a decrease in the magnitude 
of the jump). The true transition point is indicated by 
an arrow and is as far as about 25% away from the pre- 
sumed value obtained by extrapolating data to the ther- 
modynamic limit without force biases. In short, the new 
algorithm confirms that the escape transition is a nor- 
mal first-order transition with all the expected finite-size 
effects which include (approximate) the crossing of the 
curves for different size without a systematic shift and 
sharpening of the transition with increasing L. 



B. How to avoid spurious results? 



A y0 


A A A A A A A A 








L = 200, old 


new a 


L= 1600, old 


new A 



6 6.5 7 7.5 8 

N/L 

FIG. 16: Average end-to-end distance divided by L, R^/L, 
plotted against N/L for L = 200 and 600, and D = 21. Re- 
sults are obtained by PERM with (new) and without (old) 
force biases. The true transition point is indicated by an ar- 
row. 



discontinuity in the slope leaving the signature of a first- 
order transition. On the other hand, the extrapolated 
curve of R\\/L versus N/L also acquires a discontinu- 
ity in the slope suggesting a continuous transition. A 
thoughtful investigator would recall the statement about 
impossibility of phase transitions in Id systems. How- 
ever, at the transition point (Ru — L), the chain starts 
leaving the tube and claims back (at least partially) its 
three-dimensional nature! Overall, the picture of the phe- 
nomenon is very rich and thought-provoking. The only 
curve that is slightly outside the general systematic trend 
corresponds to L = 200. However, this can be clearly at- 
tributed to the stronger finite-size effect. 

Results obtained by PERM with and without force bi- 



The question posed is of a very general nature that 
any researcher has to address again and again in vari- 
ous situations. As far as simulated systems undergoing 
phase transitions are concerned, the best answer that we 
could come up with is related to the careful analysis of 
the distribution functions of the order parameter. Here 
we do not address the question of how to define the or- 
der parameter for a particular phase transition. In the 
case of the escape transition, it was shown in Sec. IV 
that the order parameter should describe the stretching 
of the confined segments of a chain, s = r/N for the 
confined state, and s = L/n for the escaped state. It 
is intuitively clear that the stretching degree is stronger, 
i.e., the value of s is larger, for a fully confined chain (in 
an imprisoned state) than that for a partially confined 
chain (in an escaped state). 

Far from the transition point the distribution func- 
tion of the order parameter P(N, L, D, s) is unimodal 
but near the transition point it has a bimodal form as 
shown in Fig. [H] and Fig. [T7] As the system size in- 
creases, i.e. L increases for a fixed vale of D = 21, the 
shape of P(N,L,D,s) changes systematically. At the 
transition point the two peaks are of equal height and 
the depth of the gap between them increases with L. 
For L = 200 the bimodal behavior of the distribution 
function P(N, L, D, s) is developed by both algorithms 
(Fig [TTa). A more close inspection shows that using 
PERM without force biases there is a cut-off at large 
s although the maxima is traced very confidently. The 
cut-off appears earlier as L increases and the problem 
of developing a right-hand maximum becomes progres- 
sively severe. Already for L = 800 the branch of the 
distribution function corresponding to the escaped state 
is completely cut off. For L = 1600 even the full descrip- 
tion of the confined state is lost, although the vicinity of 
the maximum is reproduced very fairly (here the differ- 
ence between the old and the new method is only due to 
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FIG. 17: Histograms of the order parameter P(N, L, D, s) vs. s near the transition point (JV/L)tr for D = 21 and various values 
of L. The values of N/L are (a) 6.9, (b) 6.8, (c) 6.7, and (d) 6.7. The distribution is normalized by choosing ^ s P(N, L, D, s) = 
1. The left and right branches correspond to chains in the imprisoned state and escaped states respectively. 



normalization). A deficient sampling leads to spurious 
behavior of the equilibrium average. Since the problem 
is clearly related to poor samplings of strongly stretched 
states, introducing the force biased strategy that enriches 
stretched conformations is a natural improvement of the 
algorithm PERM. 



Even if the correct distribution function was not avail- 
able, the systematic changes in the shape of the distribu- 
tion function obtained without force biases ring the alarm 
bell: the initially smooth bimodal distribution function 
is being progressively and forcibly cut off which does not 
correspond to any reasonable physical picture. This is 
a very strong indication that the algorithm experiences 
severe difficulties in sampling the relevant portions of the 
phase space. Once these portions are identified, a reason- 
able recipe for enriching the relevant set of configurations 
could be normally introduced. To finalize, we stress that 
following the behavior of the distribution functions with 
the increase in the system size is a powerful test of the 
quality of results. 



VII. SUMMARY 

In this paper, we have shown that with the force-biased 
PERM we are able to produce sufficient samplings for 
the configurations in the escaped regime for the 3d es- 
cape problem. It solves the difficulty of getting inhomo- 
geneous flower-like conformations mentioned in the 2d 
escape transition [l^|. As predicted by the Landau the- 
ory approach, we indeed see rather pronounced jumps in 
the quantities of the end-to-end distance, the fraction of 
imprisoned number of monomers, and the order param- 
eter. These jumps become sharper as the tube length L 
increases or the tube diameter D decreases, indicating 
the transition is a first-order like. The occurrence of an 
abrupt change of the slope of the free energy gives the 
precise estimate of the transition point for the polymer 
chain of length N confined in a finite tube of length L 
and diameter D. We are also able to give an evidence 
for the two minimum picture of the first-order like tran- 
sition, and our numerical results are in perfect agreement 
with the theoretical predictions by Landau theory (27| . 

For fully confined single polymer chains in a tube, 
we have also shown nice cross-over data collapse for the 
free energy and the end-to-end distance with high accu- 
racy MC data by using the algorithm PERM with fc-step 
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Markovian anticipation. Although the scaling behaviour 
is well know for this problem, it has never been shown 
directly in the literatures apart from some results already 
given in Ref. [27| . 

We expect that the method developed in the present 
paper, which allows to obtain equilibrium properties of 
chains which are partially in free space and partially con- 
fined in cylindrical tubes, could also be useful to clar- 
ify some aspects of the problem of (forced) transloca- 
tion of long polymers through narrow pores in mem- 
branes. Of course, then one needs to consider configu- 
rations of chains which have (in general) escaped parts 
on both sides of the confining tube. In this problem, it is 
strongly debated to what extent these transient, partly 
escaped, configurations requirement thermal equilibrium 
states [23, [29|, |30, |33, [32|, [33, yj, [33, [301 • Of course, for ap- 
plications to biomolecules such as RNA, single-stranded 



DNA more realistic models than self-avoiding walks on 
a lattice need to be used, before one may makes contact 
with the experiment [33, [38[ . 



Acknowledgements 

We are grateful to the Deutsche Forschungsgemein- 
schaft (DFG) for financial support: A. M.S. and L.I.K. 
were supported under grant Nos. 436 RUS 113/863/0, 
and H.-P.H. was supported under grant NO SFB 
625/A3. A. M.S. received partial support under grant 
NWO-RFBR 047.011.2005.009, and RFBR 08-03-00402- 
a. Stimulating discussions with A. Grosberg, J.-U. Som- 
mer, and Z. Usatenko are acknowledged. H.-P.H. thanks 
P. Grassberger for very helpful discussions. 



[1] G. Subramanian, D. R. M. Williams, and P. A. Pincus, 
Europhys. Lett., 29, 285 (1995); Macromolecules, 29, 
4045 (1996). 

[2] J. Ennis, E. M. Sevick, and D. R. M. Williams, Phys. 

Rev. E, 60, 6906 (1999). 
[3] E. M. Sevick and D. R. M. Williams, Macromolecules, 

32, 6841 (1999). 
[4] B. M. Steels, F. A. M. Leermakers, and C. A. Haynes, J. 

Chrom. B., 743, 31 (2000). 
[5] A. Milchev, V. Yamakov, K. Binder, Phys. Chem. Chem. 

Phys., 1, 2083 (1999); Europhys. Lett., 47, 675 (1999). 
[6] P. Grassberger, Phys. Rev. E 56, 3682 (1997). 
[7] C. M. Schroeder, H. P. Babcock, E. S. G. Shaqfeh, and 

S. Chu, Science 301, 1515 (2003). 
[8] T. Hugel, M. Seitz, Macromol. Rapid Commun. 22, 989 

(2001). 

[9] H. Matsuoka, Macromol. Rapid Commun. 22, 51 (2001). 
[10] M. P. Allen, D. J. Tildesley, "Computer Simulation of 

liquids" , Oxford University Press (1989). 
[11] A. M. Skvortsov, L. I. Klushin and F. A. M. Leermakers, 

J. Chem. Phys. 126, 024905 (2007). 
[12] L. D. Landau and M.E. Lifshitz. "Statistical Physics", 

Pergamon Press, Oxford (1958). 
[13] M. Daoud, P. G. de Gennes, J. Phys. (paris), 38, 85 

(1977). 

[14] K. Kremer and K. Binder, J. Chem. Phys., 81, 6381 
(1984). 

[15] A. Milchev, W. Paul, and K. Binder, Macromol. Theory 

Simul., 3, 305 (1994). 
[16] P. Sotta, A. Lesne, and J. M. Victor, J. Chem. Phys. 

112, 1565 (2000). 
[17] Y. Yang, T. W. Burkhardt, and G. Gompper, Phys. Rev. 

E, 76, 011804 (2007). 
[18] H.-P. Hsu, K. Binder, L. I. Klushin, and A. M. Skvortsov, 

Phys. Rev. E 76, 021108 (2007). 
[19] W. Reisner, K. J. Morton, R. Riehn, Y. M. Wang, Z. Yu, 

M. Rosen, J. C. Sturm, S. Y. Chou, E. Frey, and R. H. 

Austin, Phys. Rev. Lett. 94, 196101 (2005). 
[20] R. Randel, H. Loebl, and C. Matthai, Macromol. Theory 



Simul. 13, 387 (2004). 
[21] H.-P. Hsu and P. Grassberger, Eur. Phys. J. B, 36, 209 

(2003) . 

[22] P. G. de Gennes, "Scaling Concepts in Polymer Physics" 

; (Cornell University Press Ithaca, New York, 1979). 
[23] H.-P. Hsu, P. Grassberger, Macromolecules 37, 4658 

(2004) . 

[24] P. Grassberger, J. Phys. A: Math Gen. 38, 323 (2005). 
[25] T. W. Burkhardt and I. Guim, Phys. Rev. E 59, 5833 
(1999). 

[26] L. I. Klushin, A. M. Skvortsov, and F. A. M. Leermakers, 

Phys. Rev. E 69, 061101 (2004). 
[27] L. I. Klushin, A. M. Skvortsov, H.-P. Hsu, and K. Binder, 

Macromolecules 41, 5890 (2008). 
[28] W. Sung and P. J. Park, Phys. Rev. Lett. 77, 783 (1996). 
[29] M. Muthukumar, J. Chem. Phys. Ill, 10371 (1999). 
[30] D. K. Lubensky and D. R. Nelson, Biophys. J. 77, 1824 

(1999). 

[31] J. Chuang, Y. Kantor, and M. Kardar, Phys. Rev. E 65, 
011802 (2001). 

[32] Y. Kantor and M. Kardar, Phys. Rev. E 69, 021806 
(2004). 

[33] A. Milchev and K. Binder, J. Chem. Phys. 121, 6042 
(2004). 

[34] K. Luo, T. Ala-Nissila, and S.-C. Ying, J. Chem. Phys. 
124, 034714 (2006); K. Luo, I. Huopaniemi, T. Ala- 
Nissila and S.-C. Ying, J. Chem. Phys. 124, 114704 

(2006) . 

[35] J. K. Wolterink, G. T. Barkema, and D. Panja, Phys. 
Rev. Lett. 96, 208301 (2006); D. Panja, G. T. Barkema 
and R. C. Ball, J. Phys.: Condens. Matter 19, 432202 

(2007) ; ibid. 20, 095224 (2008). 

[36] J. L. A. Dubbeldam, A. Milchev, V. G. Rostiashvili, and 
T. A. Vilgis, Phys. Rev. E 76, 010801(R) (2007); Euro- 
phys. Lett. 79, 18002 (2007). 

[37] B. Alberts et al. Molecular Biology of the Cell (Gerland 
Publ., New York, 1994). 

[38] A. Meller, J. Phys.: Condens. Matter 15, R581 (2003). 



